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CO  S  Oi 


1.0  INTRODUCTION 

The  primary  objective  of  this  Phase  I  project  was  to  demonstrate  the  use  of  a  powerful 
technique  to  determine  hypocenter  estimates  from  the  arrival  times  of  regional  or  local  seismic 
body  waves  and  3-D  velocity  models.  We  have  focused  on  the  development  of  a  grid-search 
based,  nonlinear  hypocenter  inversion  technique  which  can  greatly  improve  estimates  of  event 
locations.  Our  algorithm  calculates  travel  times  using  a  finite  difference  approximation  of  the 
eikonal  equation  (Podvin  and  Lecomte,  1991).  The  location  algorithm  determines  the  event 
hypocenters  by  performing  a  search  over  a  3-D  grid  of  potential  hypocenters  to  minimize 
arrival  time  residuals.  We  have  modified  the  travel  time  prediction  algorithm  to  begin 
exploiting  secondary  regional  phase  travel  times  and  have  begun  optimizing  the  grid-search 
location  algorithm  for  the  regional  monitoring  problem.  In  this  paper  we  will  report  on  the 
method’s  capabilities  in  a  regional  seismic  monitoring  application  using  synthetic  data  and 
a  preliminary  3-D  velocity  model  we  have  built  for  the  Pakistan/India  region  under  other 
work  funded  by  the  Defense  Threat  Reduction  Agency. 

1.1  Background 

Current  national  monitoring  goals  require  the  accurate  location  and  characterization  of  seis¬ 
mic  events  with  a  high  degree  of  confidence.  The  global  extent  of  the  monitoring  efforts, 
which  requires  the  swift  detection  and  identification  of  large  numbers  of  small  magnitude 
events,  requires  the  development  of  new  capabilities  to  analyze  regional  seismic  data.  A  crit¬ 
ical  step  in  correctly  identifying  a  seismic  event  is  the  determination  of  an  accurate  source 
location.  Certain  regions  of  the  world  have  been  adequately  instrumented  and  calibrated 
to  produce  location  estimates  with  high  accuracy.  But  in  many  other  regions,  systematic 
errors  due  to  a  lack  of  regional  location  calibration  information  can  cause  significant  devia¬ 
tions  between  the  calculated  and  true  source  locations.  There  are  two  predominant  ways  to 
improve  the  location  of  a  seismic  event.  One  is  to  use  reliable  path  calibrations  developed 
from  events  of  known  location.  The  other  method  is  to  use  high-resolution,  3-D  models  of 
the  Earth’s  velocity  structure  for  predicting  travel  times  in  the  location  process.  The  latter 
method  is  ultimately  preferable  because  of  the  obvious  side  benefits  that  accompany  it,  such 
as  the  ability  to  use  the  3-D  model  in  characterizing  the  uncertainty  of  the  final  location 
estimate.  Models  of  the  Earth  can  never  account  for  all  the  complexity  seen  in  seismic  data, 
but  sufficiently  detailed  models  of  specific  regions  of  interest  can  be  constructed  to  produce 
locations  accurate  enough  for  regional  seismic  monitoring. 

1.2  Nonlinear  Event  Location  Using  3-D  Velocity  Models 

Finding  the  precise  location  of  an  earthquake  or  a  man-made  seismic  event  has  been  a  prob¬ 
lem  of  long-standing  interest  to  seismologists.  All  methods  of  seismic  event  location  require 
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the  prediction  of  travel  times  between  an  event  hypocenter  and  the  recording  stations,  using 
some  form  of  a  velocity  model  between  the  event  and  the  station.  Hypocentral  locations  are 
then  estimated  through  the  minimization  of  the  travel  time  residuals  between  the  observed 
and  calculated  arrival  times.  The  current  location  procedure  at  the  Prototype  International 
Data  Center  (PIDC)  uses  a  conventional  least-squares  algorithm  to  locate  seismic  events. 
In  addition  to  the  first-arriving  and  S  phases,  the  PIDC  incorporates  a  range  of  secondary 
phases  (P9,  pP,  sP,  PcP,  PP,  etc.)  as  well  as  azimuth  and  slowness  measurements  obtained 
from  both  arrays  and  three-component  stations.  However,  the  PIDC  is  currently  limited 
to  using  global  travel  time  tables  (IASP91)  in  most  regions  of  the  world  (GSE  Conference 
Room  Paper  243, 1995),  although  methods  to  incorporate  station  corrections,  regional  travel 
time  models,  and  path-dependent  measurements  are  being  tested.  Other  methods  of  seis¬ 
mic  event  location  estimation  have  been  developed  in  recent  years  which  take  advantage 
of  computational  advances  such  as  grid-based  ray  tracing  methods  (Vidale,  1988;  Moser, 
1991;  Podvin  and  Lecomte,  1991).  These  ray  tracing  methods,  combined  with  sophisticated 
inversion  techniques,  have  led  to  grid-search  earthquake  location  algorithms  capable  of  pre¬ 
dicting  hypocenters  in  3-D  media  (Moser  et  al.  1992;  Wittlinger  et  al.,  1993;  Lomax  et  al., 
1998).  The  benefits  of  these  methods  are  that  they  do  not  require  the  calculation  of  travel 
time  derivatives  near  an  estimated  hypocenter  (as  do  the  traditional  methods)  and  that  they 
are  less  susceptible  to  instabilities  in  the  numerical  solution.  Some  of  them  are  also  able 
to  handle  complicated  3-D,  regional  velocity  models,  which  are  much  more  appropriate  for 
the  modern  seismic  monitoring  problem  in  many  cases.  However,  grid-based  location  algo¬ 
rithms  still  require  improvements  before  they  can  become  routinely  useful  in  the  regional 
monitoring  problem.  During  Phase  I  of  this  research,  we  have  developed  a  version  of  a  grid- 
based  location  algorithm  tailored  to  the  regional  monitoring  setting.  We  have  improved  the 
accuracy  of  the  travel  time  prediction  algorithm  using  a  multi-grid  approach,  incorporated 
some  secondary  phases  in  the  location  estimate,  and  increased  the  efficiency  of  the  3-D  grid 
location  technique.  We  have  applied  the  new  method  to  a  synthetic  data  sets  obtained  by 
ray  tracing  through  a  preliminary  3-D  model  of  the  Pakistan  region  constructed  for  use  in 
other  work  funded  by  the  Defense  Threat  Reduction  Agency.  Each  of  these  improvements 
will  be  discussed  in  the  following  sections. 

2.0  RESEARCH  ACCOMPLISHED 

2.1  Task  1:  Finite  Difference  Computation  of  Regional  Travel  Times 

An  essential  component  of  any  location  algorithm  is  the  ability  to  predict  the  travel  times  of 
seismic  phases  through  a  given  velocity  model.  In  the  regional  monitoring  problem,  lateral 
and  3-D  heterogeneity  in  the  crust  precludes  the  use  of  traditional,  global  1-D  velocity  models 
such  as  IASP91  or  J-B.  3-D  velocity  models  require  more  sophisticated  travel  time  modeling 
routines;  thus,  we  use  a  3-D  eikonal  equation  solver  originally  developed  by  Podvin  and 
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Lecomte  (1991).  They  extended  the  original  method  of  Vidale  (1988)  in  various  ways.  The 
Podvin-Lecomte  algorithm  (hereafter  referred  to  as  P-L)  correctly  takes  into  account  the 
possibility  that  several  locally  independent  wavefronts  may  contribute  to  the  wavefield  at 
any  point  in  the  medium.  Multiple  arrivals  at  any  point  are  systematically  considered,  and  a 
first  arrival  criterion  is  used  to  pick  the  one  with  the  minimum  travel  time.  The  traditional 
eikonal  equation  methods  do  not  account  for  multi-pathing,  using  instead  a  single-wavefront 
approximation  in  the  propagation.  This  is  the  main  reason  why  the  traditional,  Vidale-style 
finite  difference  algorithm  cannot  correctly  handle  large  velocity  contrasts.  The  P-L  method 
can  handle  velocity  contrasts  as  high  as  1:10,  no  matter  what  the  shape  of  the  feature. 

In  the  P-L  algorithm,  the  velocity  model  is  set  up  as  a  uniform  Cartesian  grid  in  three 
dimensions.  The  first  arrival  travel  time  is  computed  at  each  node  of  the  grid  from  a  source 
located  within  the  grid.  In  our  application,  the  “source”  is  actually  taken  to  be  a  station  and, 
owing  to  reciprocity,  each  node  of  the  3-D  grid  is  the  hypocenter  of  some  event.  The  resulting 
grids,  one  for  each  station  and  phase,  embody  a  set  of  3-D  travel  time  tables  which  can  be 
used  by  an  event  location  algorithm.  The  ability  to  compute  travel  times  from  one  point  (i.e. 
the  station)  to  every  other  grid  point  is  a  great  advantage  that  the  eikonal  solvers  (and  other 
grid  methods)  have  over  traditional  bending  and  shooting  methods.  Travel  times  for  S  phases 
are  calculated  using  a  shear  wave  velocity  grid,  or  S  arrival  times  can  be  predicted  from  a 
pre-set  Vp/Vs  ratio.  This  method  of  travel  time  computation  is  considerably  faster  than 
two-point  ray  tracing  when  a  large  number  of  sources  or  receivers  is  involved,  distributed 
throughout  the  model. 

During  Phase  I,  we  adapted  the  P-L  algorithm  to  the  regional  seismic  monitoring  setting. 
Several  changes  in  the  technique  had  to  be  implemented  before  the  algorithm  was  suitable 
for  our  use.  A  major  concern  to  address  was  the  accuracy  of  the  method  when  used  in 
a  regional  monitoring  scenario,  where  source-receiver  distances  are  large.  For  example,  a 
1%  error  in  the  travel  time  of  P9  at  an  epicentral  1000  km  amounts  to  2  seconds  of  error, 
sufficient  to  cause  an  unacceptable  bias  in  the  location  of  an  event.  To  deal  with  this  concern, 
we  made  changes  in  the  initialization  procedure  to  greatly  increase  the  accuracy  of  the  travel 
time  estimates  in  the  immediate  region  surrounding  the  source  (i.e.  the  station,  in  our  case). 
Local  computations  in  the  P-L  method  are  based  on  a  plane  wavefront  approximation,  which 
near  the  source  can  be  quite  imprecise.  To  improve  the  accuracy  of  the  method  in  the  case 
where  the  source  is  in  a  heterogeneous  zone,  we  have  implemented  a  multi-grid  approach. 
This  entails  extracting  a  cubic  region  surrounding  the  source  and  calculating  travel  times 
on  a  grid  with  much  finer  mesh  spacing.  The  travel  times  from  the  nodes  of  the  finer  mesh 
common  to  the  coarser  grid  are  reinserted,  and  the  algorithm  proceeds  to  propagate  the 
travel  time  grid  to  all  nodes  using  finite  differences.  A  schematic  of  this  process  is  shown  in 
Figure  1. 

It  was  also  necessary  to  adapt  the  P-L  method  to  the  regional  setting  by  performing 
transformations  of  the  velocity  model  and  output  times  between  geographic  and  Cartesian 
coordinates.  The  P-L  code  was  originally  written  for  local  earthquake  location  studies,  where 
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Figure  1:  Schematic  of  the  multi-grid  process  to  increase  the  accuracy  of  the  Podvin-Lecomte 
travel  time  prediction  method. 
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curvature  of  the  Earth  could  be  safely  ignored.  Rather  than  completely  rewriting  the  P-L 
code  during  Phase  I,  we  have  implemented  flattening  transformations  on  the  velocity  model 
prior  to  running  the  code  and  then  ’’unflatten”  the  resulting  travel  time  grids  during  the 
location  procedure.  The  flattening  equations  are 
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where  /  and  v'  are  the  transformed  depths  and  velocities  that  the  P-L  code  uses.  In  addition, 
v,  r  and  z  =  a  —  r  are  the  velocities  and  depths  in  the  spherical  earth,  with  a  —  6371 km. 
This  procedure  is  exact  for  1-D  earth  models  and  approximate  for  3-D  models. 

After  making  these  changes  to  the  P-L  algorithm,  we  performed  extensive  testing  of  the 
method  to  determine  the  accuracy  for  various  velocity  models.  We  did  this  prior  to  using 
the  newly  improved  method  on  our  strongly  heterogeneous,  preliminary  3-D  model  of  the 
Pakistan  region.  The  models  we  tested  were  a  homogeneous  half-space,  a  2-layer  model 
mimicking  the  velocity  contrast  between  the  crust  and  mantle,  and  the  IASPEI91  1-D  global 
velocity  model  to  a  depth  of  410  km.  Our  goal  was  to  produce  travel  time  errors  less  than  0.5 
seconds  for  all  azimuths,  source  depths  and  epicentral  distances  to  2000  km.  The  tests  were 
performed  to  determine  algorithmic  parameters  such  as  optimal  grid  spacing  in  the  coarse 
grid,  the  overall  dimensions  of  the  source-region  initialization  grid,  and  grid  spacing  in  the 
source-region  grid.  Based  on  these  tests,  we  are  currently  using  5  km  ( h  =  5)  spacing  in  the 
coarse  regional  grid,  a  overall  source  region  size  of  22 h  per  side  with  the  source  point  at  the 
center  node,  and  a  source  region  grid  spacing  of  1  km.  In  all  of  the  tests  discussed  below, 
we  have  designed  the  models  to  support  the  work  we  are  doing  on  our  regional  Pakistan 
model.  Therefore,  we  placed  the  source  (station)  at  the  location  of  NIL,  the  IMS  station  at 
the  Earth’s  surface  in  Nilore,  Pakistan,  and  propagated  travel  times  on  a  grid  encompassing 
our  Pakistan  regional  model  (461  nodes  in  longitude  x  401  nodes  in  latitude  x  17  nodes  in 
depth.  Essentially,  the  plots  described  below  depict  the  travel  time  errors  from  sources  at 
a  given  depth  everywhere  in  the  particular  model  to  a  station  positioned  at  NIL  (33.65°N, 
73.2512°E).  A  summary  of  the  results  follows: 

1.  Homogeneous  model  —  This  model  was  produced  as  an  initial  test  of  the  method 
because  we  can  compare  the  results  to  the  simple  analytical  solution.  We  display  all 
the  results  in  this  report  as  differences  in  travel  time  between  the  P-L  method  and  the 
analytical  answer  for  horizontal  slices  through  the  resultant  travel  time  grid.  Figures 
2  and  3  are  slices  from  the  travel  time  grid  at  0  km  depth  and  60  km  depth.  The 
errors  are  all  0.4  seconds  and  under,  with  the  largest  errors  concentrated  along  certain 
azimuths.  Apparently,  the  travel  times  along  raypaths  passing  through  the  corners  or 
model  cells  are  most  accurate,  while  those  that  intersect  cell  faces  are  least  accurate. 
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MultiGrid  Travel  Time  Error  for  Homogenous  Model 
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Figure  2:  The  travel  time  error  for  a  homogeneous  model  (the  difference  between  the  P-L 
method  and  the  analytical  solution)  for  a  source  at  0  km  depth  on  the  depth  surface  z  0 
km  of  the  grid).  The  maximum  error  is  around  .4  seconds. 
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MultiGrid  Travel  Time  Error  for  Homogenous  Model 
(60  km  Depth,  5  km  spacing) 
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Figure  3:  The  travel  time  error  for  a  homogeneous  model  (the  difference  between  the  P-L 
method  and  the  analytical  solution)  for  a  source  at  0  km  depth  on  the  depth  surface  z  =  60 
km  of  the  grid).  The  maximum  error  is  slightly  over  .4  seconds. 
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2.  Two  layer  (crust /mantle)  model  with  flat  interface  —  This  model  consists  of  a 
35  km  thick  layer  of  6.0  km/s  over  a  halfspace  layer  of  8.0  km/s.  An  analytical  solution 
also  exists  for  this  example,  and  we  have  plotted  the  travel  time  surfaces  from  the  grid 
at  0  km  and  60  km  depth  (Figures  4  and  5).  The  narrow  black  ring  surrounding  the 
station  in  Figure  4  is  an  artifact  of  the  chosen  color  scale,  and  consists  of  negative 
numbers  very  close  to  zero.  It  is  related  to  the  crossover  distance  between  P9  and  P„. 
The  absolute  errors  do  not  exceed  approximately  0.3  seconds  at  any  point  in  the  grid. 


3.  IASP91  1-D  global  model  —  The  IASPEI91  model  is  one  of  the  currently  used 
default  models  for  location  processing  in  IDC  software  and  databases.  Therefore,  we 
tested  the  multi-grid  version  of  the  P-L  code  on  the  1-D  layered  model  using  a  similar 
set-up  as  with  the  previous  test  models.  One  exception  to  previous  models  was  that 
we  extended  our  Pakistan  regional  model  to  410  km  depth  by  blending  it  into  the 
IASP91  model  below  80  km,  which  is  the  maximum  depth  of  the  3-D  regional  Pakistan 
model.  Figures  6  and  7  depict  the  error  surfaces  in  the  travel  time  grid  at  20  km  and 
60  km  depth  for  a  source  sited  at  the  location  of  the  station  NIL,  as  in  the  above 
examples.  The  errors  are  with  respect  to  “analytic”  IASP91  travel  times  generated 
by  limearly  interpolating  the  IASP91  P  wave  table  (obtained  from  the  Center  for 
Monitoring  Research)  onto  the  same  3-D  grid  used  in  the  P-L  calculation.  The  2-D 
IASP91  table  is  sampled  only  every  one  degree  in  distance  and  variably  in  depth  (from 
5  to  50  km  depth  spacing),  so  the  times  are  not  exact.  This  results  in  approximately  a 
-1.3  second  discrepancy  in  our  P-L  results  at  20  km  depth  in  a  circle  around  the  source 
at  a  distance  corresponding  to  the  P9/Pn  cross-over  distance,  which  is  not  properly 
represented  in  the  1-degree  sampled  analytic  P  table.  At  60  km,  the  branch  of  the 
travel  time  vs.  distance  curve  corresponding  to  upgoing  rays  has  significant  concave- 
upward  curvature,  and  this  also  is  not  well  represented  in  the  tables.  Hence,  there  is 
approximately  a  -1.6  second  discrepancy  in  a  ring  under  the  source.  Ignoring  these 
discrepancies,  which  are  clipped  in  the  plots,  the  we  see  that  the  errors  in  our  P-L 
calculations  are  generally  less  than  0.5  seconds. 

Following  the  accuracy  tests  on  the  simpler  models,  the  P-L  algorithm  was  also  modified 
to  produce  travel  time  tables  for  secondary  phases  in  our  regional  Pakistan  model.  For  now 
we  have  concentrated  on  developing  travel  time  tables  for  P9,  Pn  and  true  mantle  P  phases, 
as  these  are  by  far  the  most  frequently  observed  body  wave  phases  in  seismic  event  bulletins. 
The  extension  to  S  waves  is  fairly  trivial,  and  we  plan  to  implement  this  capability  as  the 
research  continues.  In  the  following  paragraphs  we  describe  our  procedure  for  making  travel 
time  tables  for  P9,  Pn  and  P. 

We  adopt  the  convention  that  the  “P  wave”  travel  time  table  used  for  event  location 
should  contain,  at  each  grid  point,  the  travel  time  of  the  first  arriving  phase  to  that  point. 
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MultiGrid  (lOh)  Travel  Time  Error  for  2  Layer  Model 
(0  km  Depth,  5  km  spacing) 


55  60  65  70  75  »0 


Longitude  (deg) 


Figure  4:  The  travel  time  error  between  the  P-L  method  and  an  analytical  solution  for  a 
2-layer  model  (35  km  thick  layer  of  6.0  km/s  over  a  halfspace  of  8.0  km/s).  The  source  is 
placed  at  the  location  of  the  station  NIL,  and  the  depth  surface  shown  is  2  =  0  km  slice  of 
the  grid.  The  absolute  maximum  error  does  not  exceed  .35  seconds. 


9 


Time  (s) 


MultiGrid  (lOh)  Travel  Time  Error  for  2  Layer  Model 
(60  km  Depth,  5  km  spacing) 
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Figure  5:  The  travel  time  error  between  the  P-L  method  and  an  analytical  solution  for  a 
2-layer  model  (35  km  thick  layer  of  6.0  km/s  over  a  halfspace  of  8.0  km/s).  The  source  is 
placed  at  the  location  of  the  station  NIL,  and  the  depth  surface  shown  is  2  =  60  km  slice  of 
the  grid.  The  absolute  maximum  error  does  not  exceed  .35  seconds. 
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MultiGrid  (lOh)  Travel  Time  Error  for  IASPx  Model 


(20  km  Depth,  5  km  spacing) 
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Figure  6:  The  travel  time  error  between  the  P-L  method  and  travel  times  from  the  IASP91 
1-D  velocity  model.  The  source  is  placed  at  the  location  of  the  station  NIL,  and  the  depth 
surface  shown  is  2  =  20  km  slice  of  the  grid.  Errors  in  the  majority  of  the  grid  do  not  exceed 
.3  seconds;  the  larger  errors  near  the  station  are  due  to  the  inexact  times  produced  by  our 
IASP91  analytic  solution. 
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MultiGrid  (lOh)  Travel  Time  Error  for  IASPx  Model 
(60  km  Depth,  5  km  spacing) 


Longitude  (deg) 
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Figure  7:  The  travel  time  error  between  the  P-L  method  and  travel  times  from  the  IASP91 
1-D  velocity  model.  The  source  is  placed  at  the  location  of  the  station  NIL,  and  the  depth 
surface  shown  is  the  z  =  60  km  slice  of  the  grid.  Errors  in  the  majority  of  the  grid  do  not 
exceed  .4  seconds;  the  larger  errors  near  the  station  are  due  to  the  inexact  times  produced 
by  our  IASP91  analytic  solution. 
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The  first  arriving  phase  may  be  P9  (or  some  other  purely  crustal  phase),  Pn  (head  wave 
along  the  Moho  discontinuity),  or  mantle  P  (i.e.,  rays  turning  or  originating  below  the 
Moho),  depending  on  the  event  location  with  respect  to  the  station.  By  default,  the  P-L 
algorithm  computes  first-arrival  travel  times. 

To  compute  a  travel  time  table  for  P9  or  Pn,  whether  or  not  it  is  first,  we  modify  the 
velocity  models  input  to  the  P-L  in  certain  ways.  To  produce  P9  travel  times,  we  modify 
the  velocity  throughout  the  mantle  to  be  very  small,  such  that  the  first  arrival  to  a  surface 
station  from  any  event  in  the  crust  will  be  P9.  (Strictly,  it  will  be  a  crustal  phase  of  some 
sort,  which  is  our  operational  definition  of  P9  for  the  purposes  of  regional  event  location). 
The  travel  time  grid  produced  in  this  way  will  have  meaningless  travel  times  at  grid  points 
(hypocenters)  in  the  mantle.  We  edit  the  table  to  change  these  to  minus  1.0  seconds,  the 
signal  to  the  grid-search  location  code  that  the  travel  time  is  undefined. 

To  compute  Pn  times,  we  modify  the  model  differently.  At  each  horizontal  location  (i.e., 
one  column  of  the  velocity  model  grid)  the  mantle  velocity  is  made  constant  with  depth  and 
equal  to  the  P„  velocity,  i.e.,  model  cells  deeper  in  the  mantle  are  given  the  same  velocity 
as  the  shallowest  mantle  cell.  Furthermore,  the  flattening  transformation  is  not  applied  to 
the  mantle  cells,  which  would  introduce  a  velocity  gradient  with  depth.  In  this  way,  the 
first  arriving  phase  from  any  event  in  the  crust  will  be  either  P9  or  Pn;  rays  will  not  turn 
in  the  mantle.  The  resulting  travel  time  grids  are  then  edited  in  two  ways.  First,  points  in 
the  mantle  are  changed  to  -1.0.  Second,  if  the  travel  time  for  a  point  in  the  crust  is  greater 
than  or  equal  to  the  P9  time,  that  point  is  also  assigned  -1.0  seconds.  This  eliminates  P9 
times  from  the  Pn  table.  It  also  eliminates  P„  as  a  secondary  arrival  behind  P9,  but  this 
only  occurs  over  a  limited  distance  range  from  a  station  (between  critical  and  cross-over 
distances)  and,  in  any  case,  the  ability  to  observe  secondary  P„  is  difficult  in  practice. 

After  the  complete  set  of  phase  travel  time  tables  are  calculated  for  the  Pakistan  regional 
model,  three  (P,  P9,  P„)  per  station,  they  are  used  in  computing  the  objective  function  that 
is  optimized  by  our  grid-search  location  algorithm. 

2.2  Task  2:  Grid-Search  Location  Algorithm 

The  other  key  component  of  our  location  algorithm  besides  travel  time  computation  is  an 
inversion  method  to  estimate  the  location  of  an  event  from  observed  arrival  times.  Our 
algorithm  finds  maximum  likelihood  estimates  of  the  hypocenter,  origin  time,  and  data 
variance  scale  factor  of  the  data  errors.  Under  a  fairly  general  probability  model  of  the  errors, 
maximizing  the  likelihood  corresponds  to  minimizing  the  norm  of  arrival  time  residuals,  i.e. 
a  data  “misfit”  objective  function. 

The  hypocentral  parameters  of  a  seismic  event  are  a  three-dimensional  position  vector  x 
and  an  origin  time  t.  Let  d  =  (di,  d%. . . . ,  dn)  be  an  n-dimensional  vector  of  arrival  time  data 
from  a  seismic  network  picked  from  various  seismic  phases.  The  event  location  problem  may 
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be  expressed  as 


di  - 1  +  Tj(x)  +  e<,  i  =  l,...,n  (3) 

where  Tj  is  a  travel  time  function  (traveltime  table)  for  the  ith  datum  and  e*  is  an  error  with 
respect  to  this  function.  The  index  i  counts  over  both  stations  and  phase  types  (Pg,  Pn, 
etc.),  including  only  those  combinations  that  have  been  actually  observed  (i.e.,  completeness 
of  data  coverage  is  not  assumed). 

A  likelihood  function  is  determined  by  the  joint  probability  density  function  (p.d.f.)  of 
the  arrival  time  data,  which  in  turn  is  determined  by  a  probability  model  of  the  data  errors. 
Our  formulation  assumes  that  the  errors  are  statistically  independent  and,  following  Billings 
et  al.  (1994),  that  each  is  distributed  with  a  “generalized  Gaussian”  p.d.f.  For  simplicity  we 
consider  here  only  the  Gaussian  case.  Further,  following  earlier  formulations  of  the  event 
location  problem  (e.g.  Jordan  and  Sverdrup,  1981),  we  assume  that  the  error  variances,  , 
are  known  in  a  relative  sense  and  write 

(Ji  =  avi  (4) 

where  the  i/*  are  given  but  the  universal  scale  parameter,  a,  is  unknown.  Under  these 
assumptions,  the  likelihood  function  is  given  by 

£(x,t,g;d)  =  (2^g1-nrgl,ie>:p{~2^'1'(x’i;d)}  (6) 

where  is  a  ‘data  misfit’  function: 

*  (x,  i;  d)  =  £  {di  - 1  -  Tj(x))2 /{vi)2.  (6) 

i= 1 

The  maximum  likelihood  estimates  of  the  unknown  parameters  are  the  values  that  max¬ 
imize  L.  We  denote  these  estimates  as  xml,  tm\  and  a**.  The  maximization  may  be  sub¬ 
jected  to  prior  constraints  on  the  parameters.  In  our  current  implementation,  we  assume 
hard  bounds  on  focal  depth  (z)  and  the  scale  variance: 

0  <  z  <  zmax 

<7min  <  a  <  <rmax. 

We  place  no  restrictions  on  latitude,  longitude,  or  origin  time. 

Our  grid  search  algorithm  obtains  the  maximum  likelihood  estimates  of  the  hypocentral 
parameters  (xml  and  tml)  and  the  p.d.f.  scale  parameter  (aml),  as  defined  above.  The  algo¬ 
rithm  computes  a  “reduced”  likelihood  function,  L,  at  each  point  of  a  3-D  grid  of  hypocenters. 

This  function  is  defined  as  __ 

L(x;  d)  =  max  L(x,  t,  o\  d).  (7) 
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With  x  fixed,  maximizing  L  with  respect  to  origin  time  t  corresponds  to  minimizing  Vk 
with  respect  to  t,  which  is  achieved  by  setting  t  to  a  weighted  sum  of  the  arrival  residuals, 
di  —  Tj(x).  The  maximizing  a  is  a  weighted  r.m.s.  residual. 

Following  previous  workers  (e.g.  Dreger  et  al. ,  1998),  we  construct  the  search  grid  of 
hypocenters  as  a  set  of  sub-grids  of  increasing  fineness.  In  our  algorithm,  each  grid  is  in 
a  cylindrical  coordinate  system,  i.e.  sampled  in  distance  and  azimuth  from  some  point  and 
in  depth.  The  azimuth  sampling  varies  with  distance  to  approximately  equalize  radial  and 
angular  sampling  in  kilometers.  The  first  6  grids  used  span  approximately  6000  km  in 
epicentral  distance  and  are  distributed  to  cover  the  entire  globe  (with  some  overlap)  and 
focal  depths  from  0  to  750  km.  The  best  hypocenter  (largest  value  of  L )  from  these  six 
grids  is  taken  as  the  center  a  finer  cylindrical  grid  of  smaller  aperture.  After  searching  this 
finer  grid,  the  algorithm  generates  a  yet  smaller,  finer  grid  centered  on  the  best  hypocenter 
found  thus  far,  and  so  on.  The  final  maximum  likelihood  hypocenter,  xml,  is  taken  as  the 
hypocenter  achieving  the  largest  reduced  likelihood  amongst  all  the  points  of  all  the  grids 
searched.  In  the  examples  shown  below,  we  used  eight  grid  refinements  in  the  recursive 
search,  with  grid  spacing  starting  at  100  km,  laterally  and  in  depth,  and  shrinking  to  0.1 
km  spacing  for  the  final  and  densest  grid.  Since  the  examples  were  run  with  the  regional 
Pakistan  model,  the  global  grids  were  bypassed:  the  largest  grid  was  3000  km  in  radius, 
extended  well  beyond  the  border  of  our  model. 

Our  grid-search  algorithm  was  initially  written  to  accommodate  only  a  spherically  sym¬ 
metric  earth  model.  In  this  case,  the  travel  time  function,  Tj(x),  is  obtained  by  interpolating 
a  2-D  travel  time  table,  which  samples  travel  time  on  a  grid  of  epicentral  distances  and  focal 
depths.  We  typically  use  the  IASP91  tables  for  this  purpose.  For  a  3-D  earth  model,  Tf(x) 
is  evaluated  by  interpolating  a  3-D  travel  time  table,  i.e.  traveltime  sampled  on  a  3-D  grid. 
We  extended  our  grid-search  algorithm  to  use  3-D  tables  with  the  following  features: 

•  Since  the  3-D  tables  are  computed  on  a  Cartesian  grid,  with  depth  stretched  via  the 
flattening  transformation,  we  applied  the  following  steps  before  interpolating  the  travel 
time  for  a  given  event.  First,  the  grid  depths  stored  with  the  3-D  table  were  “unflat¬ 
tened”  into  true  depths.  Second,  the  relative  location  between  a  station  and  an  event, 
in  terms  of  spherical  distance  and  azimuth,  is  converted  to  an  approximate  relative 
location  in  Cartesian  coordinates.  This  relative  location  is  then  used  to  locate  the 
event  within  the  P-L  grid  for  each  station  and  phase. 

•  Given  the  event  location  in  Cartesian  coordinates,  travel  time  is  determined  by  tri- 
linear  interpolation  between  eight  nodes.  If  any  of  the  nodes  has  a  travel  time  of -1.0, 
the  travel  time  is  treated  as  undefined,  e.g.  P9  for  an  event  in  the  mantle. 

•  The  travel  time  tables  are  quite  large.  For  example,  one  P  wave  table  for  the  Pakistan 
model  even  when  resampled  laterally  at  20  km  and  truncated  to  200  km  in  depth  (as 
described  below),  contains  about  500,000  numbers.  Therefore,  our  algorithm  is  written 
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Table  1:  Stations  Used  for  Synthetic  Data  Tests 


Name 

Latitude 

Longitude 

NDI 

28.6833 

77.2167 

QUE 

30.1883 

66.9500 

DRP 

31.7440 

70.2028 

SARP 

31.9215 

72.6718 

SBDP 

32.2997 

70.8072 

THW 

32.7943 

71.7427 

NIL 

33.6500 

73.2517 

CHOP 

33.6583 

73.2638 

CEP 

33.8235 

71.9090 

MAIO 

36.3000 

59.4945 

ASH 

37.9500 

58.3500 

KAT 

39.2000 

56.2667 

KSH 

39.5167 

75.9731 

FRU 

42.8333 

74.6167 

Description 
New  Delhi  (Delhi) 

Quetta 

Drazinda 

Sargodha 

Shaikh  Budin 

Thammewali 

Nilore 

Chirah  Chowk 

Cherat 

Mashhad 

Ashgabat  (Ashkhabad) 
Gyzylarbat  (Kizyl-Arvat) 
Kashi  (Kashgar) 

Bishkek  (Frunze) 


to  store  only  a  portion  of  a  table  at  a  time  in  a  series  of  cache  memory  arrays.  The 
appropriate  table  file  is  revisited  when  a  new  portion  of  a  table  is  needed,  overwriting 
the  oldest  cache  array. 

2.3  Task  3:  Application  to  Synthetic  Data 

We  applied  our  grid-search  location  algorithm  to  various  hypothetical  events  in  our  Pakistan 
model  region.  The  objective  of  these  numerical  experiments  was  mainly  to  test  that  the 
extension  of  our  location  algorithm  to  3-D  travel  time  tables  was  correct.  The  examples  also 
serve  to  demonstrate  the  need  for  calibrated  3-D  travel  time  information  in  areas  of  complex 
crustal  structure  and  the  location  biases  that  are  induced  when  global  travel  time  tables  are 
used  in  such  areas. 

For  each  hypothetical  event,  noise-contaminated  synthetic  arrival  times  were  generated 
at  14  stations  in  the  model  region  (see  Table  1),  and  then  the  location  algorithm  was  applied 
to  these  data  to  estimate  the  “true”  location.  The  synthetic  data  were  calculated  from  the 
3-D  travel  time  tables  generated  by  running  the  P-L  algorithm  on  our  3-D  Pakistan  model. 
Each  synthetic  arrival  time  was  the  sum  of  the  model-based  travel  time,  as  interpolated  from 
the  appropriate  3-D  table,  an  assumed  origin  time,  and  a  zero-mean  pseudo-random  number 
to  simulate  a  “picking  error.” 

The  3-D  travel  time  tables  for  this  purpose  were  generated  with  a  single-grid  version  of 
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Table  2:  Location  Errors  for  Crustal  Events 


True  location 

No.  of 

Epicenter  error  (km) 

Depth  error  (km) 

lot. 

Ion. 

depth 

data 

2-D 

3-D 

2-D 

3-D 

30 

60 

15 

16 

74.0 

2.6 

-15.0 

2.9 

30 

65 

15 

19 

97.6 

6.9 

-15.0 

6.3 

30 

70 

15 

21 

24.9 

2.0 

-10.2 

0.2 

30 

75 

15 

22 

17.7 

1.3 

-15.0 

-1.0 

35 

60 

15 

16 

13.1 

2.9 

-15.0 

-2.4 

35 

65 

15 

24 

10.7 

1.9 

-15.0 

0.1 

35 

70 

15 

22 

25.5 

2.2 

-15.0 

2.4 

35 

75 

15 

21 

21.9 

1.1 

-15.0 

0.2 

40 

60 

15 

16 

13.7 

1.8 

-15.0 

-0.9 

40 

65 

15 

17 

18.6 

3.7 

-15.0 

0.1 

40 

70 

15 

20 

27.1 

1.6 

-15.0 

0.9 

40 

75 

15 

18 

28.6 

4.9 

-15.0 

8.3 

P-L  using  a  grid  spacing  of  5  km.  Two  tables  for  each  station  were  generated:  P9  and  first 
arrival  P.  In  the  raytracing,  the  grids  extended  to  a  depth  of  410  km  and  80  km  for  P  and  P 9, 
respectively.  The  3-D  tables  output  by  P-L  were  resampled  and  truncated  in  depth  before 
using  them  in  the  grid-search  algorithm.  They  were  resampled  with  a  horizontal  spacing  of 
20  km  laterally;  the  depth  spacing  of  5  km  was  retained.  The  P  tables  were  truncated  to  a 
maximum  depth  of  200  km.  The  80  km  depth  maximum  for  the  P9  tables  was  retained.  By 
truncating  the  P  tables  to  200  km,  we  eliminate  event  focal  depths  for  which  the  raypath  to 
one  or  more  stations  might  dive  below  the  410  limit  used  in  the  raytracing. 

The  synthetic  data  set  for  a  given  event  comprised  P  and  P9  times  at  a  subset  of  the 
14  stations  in  Table  1.  The  subset  was  determined  by  some  rules  stating  when  a  phase  is 
observed  as  a  function  of  source-receiver  distance  (A)  and  event  depth.  For  crustal  events, 
we  arbitrarily  assumed  that  P  was  observed  for  A  >  2.5  degrees,  and  P 9  was  observed  when 
A  <  7.5  degrees.  For  mantle  events,  P9  is  never  observed  and  P  is  observed  at  all  distances. 
The  noise  standard  deviation  was  set  to  0.5  sec  for  P  and  0.25  sec  for  P9. 

While  the  synthetic  data  were  generated  with  3-D  travel  time  tables,  we  performed  the 
location  of  each  event  two  ways:  using  the  same  3-D  tables,  and  using  the  2-D IASP91  tables. 
Table  2  summarizes  the  event  mislocations  resulting  each  way  for  events  at  a  depth  of  15 
km.  Twelve  such  events  were  assumed,  on  a  5  by  5  degree  grid  over  the  model  region.  The 
leftmost  columns  of  the  table  gives  the  “true”  hypocenter  of  each  event.  The  number  of  data 
used  in  the  location  of  the  event  is  next,  followed  by  the  location  errors.  The  columns  labeled 
“2-D”  are  the  errors  when  IASP91  tables  are  used  in  the  grid  search,  and  those  labeled  “3-D” 
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Table  3:  R.M.S.  Residuals  for  Crustal  Events 


True  Location 

No.  of 

RMS  Resid.  (s) 

lot. 

Ion. 

depth 

Data 

2-D 

3-D 

30 

60 

15 

16 

1.091 

0.299 

30 

65 

15 

19 

2.476 

0.235 

30 

70 

15 

21 

1.712 

0.256 

30 

75 

15 

22 

1.480 

0.286 

35 

60 

15 

16 

2.055 

0.261 

35 

65 

15 

24 

2.343 

0.271 

35 

70 

15 

22 

1.463 

0.235 

35 

75 

15 

21 

2.377 

0.243 

40 

60 

15 

16 

2.007 

0.291 

40 

65 

15 

17 

2.561 

0.255 

40 

70 

15 

20 

2.202 

0.274 

40 

75 

15 

18 

4.491 

0.251 

are  when  the  3-D  Pakistan  tables  are  used.  We  see  immediately  that  the  location  errors  are 
much  smaller  when  the  correct  3-D  tables  are  used.  The  horizontal  mislocation  errors  are 
then  very  small,  the  largest  being  about  6.9  km  and  most  less  than  3  km.  The  errors  in 
focal  depth  are  also  generally  small,  the  largest  being  8.3  km  and  many  less  than  one  km. 
Since  the  correct  tables  were  used  in  this  case,  the  mislocations  are  due  only  to  the  noise, 
i.e.  picking  errors,  in  the  cdata.  Given  the  size  of  the  picking  errors  and  the  large  number 
of  data  at  local  and  regional  distances,  it  is  expected  that  the  picking  errors  do  not  induce 
large  event  mislocations,  which  is  what  we  see.  However,  turning  now  to  the  columns  headed 
“2-D”,  we  see  that  using  IASP91  tables  in  the  grid-search  does  induce  large  errors,  despite 
the  generally  good  coverage.  Our  Pakistan  model  differs  from  the  IASP91  model  by  up  to 
several  seconds  in  travel  time,  which  we  see  induces  large  epicentral  errors  varying  from  10 
to  100  km  amomg  the  events.  Many  of  the  errors  well  exceed  the  desired  18  km  for  CTBT 
monitoring.  The  errors  in  focal  depth  resulting  from  using  IASP91  tables  are  uniformly 
negative,  with  all  but  one  located  at  the  top  of  the  search  grid,  i.e.  surface  focus. 

Another  comparison  of  the  2-D  and  3-D  locations  is  given  in  Table  3.  This  shows  the 
r.m.s.  arrival  time  residual  resulting  from  each  location  run.  When  the  correct  3-D  tables 
are  used,  they  vary  between  0.2  and  0.3  sec,  which  is  consistent  with  the  assumed  noise  level. 
For  the  IASP91  locations,  they  are  many  times  larger  in  every  case,  varying  from  1  to  over 
4  seconds. 

Tables  4  and  5  show  analogous  results  for  12  hypothetical  events  at  a  focal  depth  of  100 
km.  Now  each  data  set  contains  exactly  one  mantle  P  wave  at  each  of  the  14  stations.  The 
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Table  4:  Location  Errors  for  Mantle  Events 


True  location 

No.  of 

Epicenter  error 

(km) 

Depth  error  (km) 

lot. 

Ion. 

depth 

data 

2-D 

3-D 

2-D 

3-D 

30 

60 

100 

14 

18.4 

4.3 

7.1 

-5.3 

30 

65 

100 

14 

13.9 

2.9 

17.0 

2.5 

30 

70 

100 

14 

11.3 

1.9 

23.6 

-5.7 

30 

75 

100 

14 

14.8 

1.9 

19.4 

-8.7 

35 

60 

100 

14 

28.7 

2.7 

0.0 

-1.5 

35 

65 

100 

14 

21.3 

0.7 

-10.2 

-22.6 

35 

70 

100 

14 

9.9 

3.2 

4.6 

-2.4 

35 

75 

100 

14 

8.1 

3.0 

6.3 

-3.8 

40 

60 

100 

14 

12.6 

4.0 

-12.3 

-0.9 

40 

65 

100 

14 

23.0 

2.8 

0.0 

-50.2 

40 

70 

100 

14 

17.3 

1.1 

-49.5 

-90.1 

40 

75 

100 

14 

7.1 

3.5 

-100.0 

9.9 

use  of  3-D  tables  in  the  location  clearly  yields  smaller  mislocations,  but  as  dramatically  as 
for  shallow  events.  This  particularly  true  for  the  depth  mislocations,  which  for  a  few  events 
are  actually  larger  when  the  3-D  tables  are  used.  We  attribute  this  to  the  fact  that  focal 
is  poorly  constrained  in  any  case  for  deep  events  below  the  network.  As  before,  The  r.m.s. 
residuals  are  consistently  larger  2-D  tables  are  used,  but  not  by  as  large  a  factor. 


3.0  CONCLUSIONS 

Our  focus  during  Phase  I  of  this  research  project  was  to  demonstrate  the  capabilities  of 
a  powerful  seismic  event  location  method  which  uses  3-D  velocity  models  and  body  wave 
arrival  times  to  estimate  hypocenters.  Experts  in  seismic  monitoring  have  stated  that  the  de¬ 
velopment  of  methods  to  improve  hypocenter  estimation  for  the  regional  seismic  monitoring 
problem  is  a  research  priority  of  great  importance  (National  Research  Council  Committee 
on  Seismology  report,  1997).  Our  research  directly  contributes  to  that  effort  through  the 
development  of  a  sophisticated  technique  to  produce  improved  event  locations  from  regional 
travel  times  and  3-D  velocity  models.  The  Phase  I  effort  has  produced  a  prototype  algo¬ 
rithm  for  locating  regional  events  with  3-D  velocity  models.  This  type  of  location  method  is 
not  currently  the  norm  at  US  monitoring  agencies.  The  method  incorporates  sophisticated, 
grid-based  ray  tracing  techniques  and  search  algorithms  to  estimate  event  hypocenters.  We 
have  made  significant  progress  on  several  technical  objectives:  1)  Improvement  in  the  accu¬ 
racy  of  the  3-D,  Podvin-Lecomte  grid-based  travel  time  code  and  adaptation  to  the  regional 
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Table  5:  R.M.S.  Residuals  for  Mantle  Events 


True  Location 

No.  of 

RMS  Resid.  (s) 

lat. 

Ion. 

depth 

Data 

2-D 

3-D 

30 

60 

100 

14 

1.090 

0.285 

30 

65 

100 

14 

0.782 

0.314 

30 

70 

100 

14 

0.720 

0.305 

30 

75 

100 

14 

0.776 

0.299 

35 

60 

100 

14 

1.197 

0.316 

35 

65 

100 

14 

0.954 

0.275 

35 

70 

100 

14 

0.655 

0.270 

35 

75 

100 

14 

0.681 

0.266 

40 

60 

100 

14 

1.688 

0.321 

40 

65 

100 

14 

1.124 

0.272 

40 

70 

100 

14 

0.747 

0.319 

40 

75 

100 

14 

0.654 

0.265 

monitoring  problem;  2)  Adaptation  of  the  maximum  likelihood  grid-based  location  method 
to  the  3-D  regional  location  problem;  3)  Incorporation  of  secondary  phases  in  hypocenter 
estimates;  4)  Application  to  a  set  of  synthetic  events  from  a  region  of  US  monitoring  concern 
to  demonstrate  the  feasibility  and  potential  of  the  method.  There  are  still  significant  im¬ 
provements  and  obstacles  to  be  addressed  before  the  technique  can  be  used  in  an  operational 
setting,  and  we  will  address  these  in  our  forthcoming  Phase  II  proposal.  However,  we  believe 
that  the  technique  has  significant  potential  to  be  useful  in  both  monitoring  seismology  and 
passive  underground  facility  characterization. 
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